
## @knitr V-DEM_EXTRA_SETUP

rm(list=ls())


# from \url{https://www.r-bloggers.com/identifying-the-os-from-r/}
get_os <- function(){
  sysinf <- Sys.info()
  if (!is.null(sysinf)){
    os <- sysinf['sysname']
    if (os == 'Darwin')
      os <- "osx"
  } else { ## mystery machine
    os <- .Platform$OS.type
    if (grepl("^darwin", R.version$os))
      os <- "osx"
    if (grepl("linux-gnu", R.version$os))
      os <- "linux"
  }
  tolower(os)
}

get_machine  <- function(){
    sysinf  <- Sys.info()
    if (!is.null(sysinf)){
        my.machine  <- sysinf['nodename']
    }
    tolower(my.machine)
}


if(get_os()=="windows" & get_machine()=="mompracen") {
    path <- setwd('C:/Users/giaco/Documents/MEGAsync/')
} else if (get_os()=="windows" & get_machine()!="mompracen") {
    path <- setwd('C:/Users/grieco/Dropbox/')
} else if (get_machine()=="aramis") {
    path <- setwd('/home/giacomo68/Dropbox/')
} else {
    path <- setwd('/home/giacomo/Dropbox/')
}


library(survival)
library(xtable)
library(ggplot2)
library(english)
library(Publish)
library(showtext)
library(extrafont)
library(mice)
library(mitools)
library(sandwich)
library(lmtest)
library(magick)

Data <- read.csv('gc-jg-project/kampala/data/preliminary-data-kampala.csv',
                 stringsAsFactors=FALSE)
Data$v2x_egaldem.lag <- Data$v2x_egaldem.lag*10


Data$country_id <- as.numeric(factor(Data$country_name))

plot1 <- image_read("gc-jg-project/kampala/paper/figures/polyarchy-plot.png")
plot2 <- image_read("gc-jg-project/kampala/paper/figures/partipdem-plot.png")
plot3 <- image_read("gc-jg-project/kampala/paper/figures/delibdem-plot.png")
plot4 <- image_read("gc-jg-project/kampala/paper/figures/egaldem-plot.png")


combined_image <- image_append(
  c(
    image_append(c(plot1, plot2), stack = FALSE),  # Row 1
    image_append(c(plot3, plot4), stack = FALSE)   # Row 2
  ),
  stack = TRUE
)


image_write(combined_image,"gc-jg-project/kampala/paper/figures/app-figure-a5.png")


m <- 10

imp <- mice(Data[,c("kampala",
                 "sidea.1991",
                 "sideb.1991",
                 "my.mil.cap.lag",
                 "v2x_libdem.lag",
                 "gdp.cap.lag",
                 "pop.lag",
                 "common.law",
                 "year",
                 "rule.law.lag",
                 "aid.recipient",
                 "aid.recipient.q",
                 "All_NGO.lag",
                 "polity2.lag",
                 "country_id")], m = m, maxit = 10, seed = 3210)

models <- list()
model_results <- list()
merged_imputations <- list()
m1.model_results <- list()
m2.model_results <- list()
m3.model_results <- list()
m4.model_results <- list()
m5.model_results <- list()
m6.model_results <- list()

for (i in 1:m) {
    complete_data <- complete(imp, i)

    complete_data$.imp <- i
    
    merged_data <- merge(complete_data, 
                         Data[,c("country_id","year","country_name",
                                 "kampala.t0","kampala.t1",
                                 "sidea.1945",
                                 "sideb.1945",
                                 "sidea.force.1945",
                                 "sideb.force.1945",
                                 "sidea.no.force.1945",
                                 "sideb.no.force.1945",
                                 "sidea.force.1991",
                                 "sideb.force.1991",
                                 "sidea.no.force.1991",
                                 "sideb.no.force.1991"
                                 )], 
                        by = c("country_id", "year"),  # Note: quoted variable names
                        all.x = TRUE)  # Keep all rows from complete_data
    
    merged_imputations[[i]] <- merged_data
        

    current_data <- merged_imputations[[i]]

    m1 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
              sidea.1991
            + sideb.1991
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_id)
           ,data=current_data)

     m1.vcov_cluster <- vcovCL(m1, cluster = current_data$country_id)

     m1.model_results[[i]] <- list(
        coefficients = coef(m1),
        variance = m1.vcov_cluster,
        n = nrow(current_data)
       )
    
}

t.n.obs <- m1$n
t.n.fail <- m1$nevent
t.n.fail.tot <- table(Data$kampala)[2]
t.miss.vals <- dim(Data)[1]-t.n.obs
